{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sisl\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First TranSiesta bias example.\n",
    "\n",
    "In this example, we will take the system from [TS 1](../TS_01/run.ipynb) and perform bias calculations on it. Note, however, that applying a bias to a *pristine* bulk system is non-physical and we recommend *NEVER* doing this. TranSiesta will *not* warn you about this and will happily calculate the non-equilibrium density for any bulk system with an applied bias. This is an ***extremely*** important point and once complete with this example, you should carefully think through why this is the case.\n",
    "\n",
    "Bias calculations are very heavy because of a full DFT+NEGF calculation *per bias-point*.\n",
    "\n",
    "We will begin with creating the structures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "graphene = sisl.geom.graphene(1.44, orthogonal=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elec = graphene.tile(2, axis=0)\n",
    "elec.write('STRUCT_ELEC.fdf')\n",
    "device = elec.tile(3, axis=0)\n",
    "device.write('STRUCT_DEVICE.fdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "In this exercise, you will be familiarized with the input options that define a bias calculation. The input options are *extremely elaborate*, yet they require little intervention when using default parameters.  \n",
    "As this is your first example of running TranSiesta with applied bias there are a few things you should know:\n",
    "\n",
    "1. Do not start by performing any $V\\neq0$ calculations until you are perfectly sure that your $V=0$ calculation is well converged and well-behaved, i.e. a small `dQ` (see [TS 1](../TS_01/run.ipynb)).\n",
    "2. When performing bias calculations, you are recommended to create a new directory for each bias: `TS_<V>`.\n",
    "3. *Any* bias calculation should be a restart by using the ***closest*** bias calculation TranSiesta density matrix. This can be ensured by copying the `siesta.TSDE` file from the ***closest*** bias calculation to the current simulation directory. I.e.\n",
    "\n",
    "   - First run $V=0$ in `TS_0`, ensure convergence etc.\n",
    "   - Second run $V=0.5\\,\\mathrm{eV}$ in `TS_0.5`, copy `TS_0/siesta.TSDE` to `TS_0.5/` and start run.\n",
    "   - Third run $V=1\\,\\mathrm{eV}$ in `TS_0.5`, copy `TS_0.5/siesta.TSDE` to `TS_1.0/` and start run.\n",
    "   - etc.\n",
    "\n",
    "   For instance, when doing negative bias points, you would copy from $0$ to $-0.5\\,\\mathrm{eV}$ and so on (this example does not calculate the negative bias-points).\n",
    "\n",
    "4. All the commands required for this example can be executed like this:\n",
    "\n",
    "    ```\n",
    "    siesta --electrode RUN_ELEC.fdf > ELEC.out\n",
    "    cd TS_0\n",
    "    cp ../C.psf .\n",
    "    siesta ../RUN.fdf > TS.out\n",
    "    # Check that the charge is converged etc.\n",
    "    cp siesta.TSDE ../TS_0.5/\n",
    "    cd ../TS_0.5\n",
    "    cp ../C.psf .\n",
    "    siesta -V 0.5:eV ../RUN.fdf > TS.out\n",
    "    # Check that it has converged...\n",
    "    cp siesta.TSDE ../TS_1.0/\n",
    "    cd ../TS_1.0\n",
    "    cp ../C.psf .\n",
    "    siesta -V 1:eV ../RUN.fdf > TS.out\n",
    "    # Check that it has converged...\n",
    "    ```\n",
    "  After every calculation, go through the output to ensure everything is well-behaved. Note that the output of a bias calculation is different from a non-bias calculation, they are more detailed.\n",
    "\n",
    "5. An additional analysis (before going to the transport calculations) is to calculate the potential drop in the junction. In sisl this is easy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v0 = sisl.Grid.read('TS_0/ElectrostaticPotential.grid.nc')\n",
    "vd = (sisl.Grid.read('TS_0.5/ElectrostaticPotential.grid.nc') - v0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`vd` then contains the potential profile (in eV). To save it as a linear average bias file (remember transport is along the first lattice vector) you can execute the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vd = vd.average(1).average(2)\n",
    "dv = (vd.dcell[0, :] ** 2).sum() ** .5\n",
    "dx = dv * np.arange(vd.shape[0])\n",
    "sisl.io.tableSile('potential_0.5.dat', 'w').write_data(dx, vd.grid[:, 0, 0])\n",
    "# You can plot this using:\n",
    "#plt.plot(dx, vd.grid.ravel())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This completes all non-equilibrium calculations for this example. However, we have only calculated the non-equilibrium density and, thus, the non-equilibrium Hamiltonian. We still need to calculate the transport properties for all bias'. Basically we can only calculate the transport properties at the calculated bias values, but generally we are interested in a full $I(V)$ curve.  \n",
    "As a user, one has three options:\n",
    "\n",
    "1. Calculate $I(V)$ for the calculated biases $V$ and perform an interpolation of $I(V)$, or\n",
    "2. Interpolate the Hamiltonian to calculate $I(V)$ for all the required biases, or\n",
    "3. Calculate all non-equilibrium Hamiltonians.\n",
    "\n",
    "The first option is by far the fastest and easiest with a sometimes poor accuracy; the second option is relatively fast, and drastically improves the accuracy; while the last option is the most accurate but may sometimes be non-feasible due to insufficient computational resources.\n",
    "\n",
    "In the following, we will calculate all transmissions using option 2. Look in the manual for the options regarding the interpolation (there are two interpolation methods).  \n",
    "Go through `RUN.fdf` and find the respective block that tells TBtrans to interpolate the Hamiltonian, also notice how the energy-grid is defined in TBtrans. You will notice that this is the fastest way to calculate the $I(V)$ curve for *any* bias, it however, will not calculate any physical quantities outside the bias window.\n",
    "\n",
    "Now complete the exercise by running TBtrans for $V\\in\\{0, 0.1,  \\dots, 1\\}$ eV. Note that instead of changing the applied bias in the fdf-file, one can do:\n",
    "\n",
    "    tbtrans -V 0.4:eV RUN.fdf\n",
    "    \n",
    "to apply $V=0.4\\,\\mathrm{eV}$, *any* fdf-flag specified on the command line has precedence! The `:` is to denote an effective space, otherwise you will have to encapsulate in quotation marks `tbtrans -V \"0.4 eV\" RUN.fdf`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do not want to run the commands manually, you may use this loop command:\n",
    "```\n",
    "for V in $(seq 0 0.1 1) ; do\n",
    "   d=TBT_${V//,/.}\n",
    "   mkdir $d\n",
    "   cd $d\n",
    "   tbtrans -V \"${V//,/.}:eV\" ../RUN.fdf > TBT.out\n",
    "   cd ../\n",
    "done\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TIME**: A remark on this exercise. Think why applying a bias to a bulk system is wrong. If you can't immediately figure this out, try and create a longer system by replacing `device = elec.tile(3, axis=0)` with, say: `device = elec.tile(6, axis=0)` and redo the calculation for a given bias. Then compare the potential profiles. Furthermore, try and imagine a physical device with the potential profile you create? Would it be possible?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the transmissions\n",
    "\n",
    "Calculate the current for all $V$, then plot it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.arange(0, 1.05, 0.1)\n",
    "I = np.empty([len(V)])\n",
    "for i, v in enumerate(V):\n",
    "    I[i] = sisl.get_sile('TBT_{:.1f}/siesta.TBT.nc'.format(v)).current()\n",
    "plt.plot(V, I * 1e6); \n",
    "plt.xlabel('Bias [V]'); plt.ylabel(r'Current [$\\mu\\mathrm{A}$]');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Why is the current $0$ for $V<0.8$?  \n",
    "*Hint*: see [TB 1](../TB_01/run.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
